IntData <- function(object, data, x, intvar, intval, valname,
  continuous = "mean", discrete = "mode", n = 100L, FUN = "glm", ...) {

  if (inherits(object, "formula") && !missing(data)) {
    FUN <- match.fun(FUN)
    m <- FUN(formula = object, data = data, ...)
  ## likely not the best way to tell if there is a predict method
  } else if (is.data.frame(object$model) &&
    is.function(getS3method("predict", class(object)))) {
    m <- object
  } else stop("'object' must be a formula or an object with an appropriate
    method for 'predict' and a copy of the data")

  mdata <- m$model
  names <- all.vars(as.formula(m))

  vals <- NormDF(mdata, continuous = continuous, discrete = discrete)

  colnames(vals) <- names
  y <- names[1L]

  if (missing(intval)) {
    b <- c(mean(mdata[, intvar], na.rm = TRUE),
      sd(mdata[, intvar], na.rm = TRUE))
    intval <- rep(b[1L] + c(b[2L], 0, -b[2L]), each = n)
    valname <- c("+1 SD", "Mean", "-1 SD")
  } else intval <- rep(intval, each = n)

  newdat <- do.call(rbind, rep(list(vals), length(intval)))
  newdat[, intvar] <- intval
  minx <- min(mdata[, x], na.rm = TRUE)
  maxx <- max(mdata[, x], na.rm = TRUE)
  rx <- seq(from = minx, to = maxx, length.out = n)

  newdat[, x] <- rep(rx, length(intval)/n)

  preddat <- data.frame(predict(m, newdat, type = "response"),
    newdat[, x], factor(rep(valname, each = n), levels = valname))
  colnames(preddat) <- c(y, x, intvar)

  SS <- SimpleSlopes(object = m, x = x, intvar = intvar,
    intval = intval, values = vals)

  output <- list(ModelData = mdata, PredictedData = preddat, SS = SS)
  class(output) <- "IntData"
  attr(output, "Variables") <- c("x" = x, "y" = y, "intvar" = intvar)
  return(output)
}

## tmp <- IntData(lm(mpg ~ hp * wt, data = mtcars), x = "hp", intvar = "wt")
## plot(tmp)
## tmp[c("coefficients", "model")]

plot.IntData <- function(x, y, ...) {
  name <- attr(x, "Variables")
  p <- ggplot(x[["ModelData"]], aes_string(x = name["x"], y = name["y"])) +
    geom_line(aes_string(x = name["x"], y = name["y"],
      linetype = name["intvar"]), data = x[["PredictedData"]], ...)
  text
  return(p)
}

coef.IntData <- function(object, ...) {
  if (!is.null(object[["SS"]])) {
    object[["SS"]]["coefficients"]
  } else cat("No data to extract coefficients from", fill = TRUE)
}

vcov.IntData <- function(object, ...) {
  if (!is.null(object[["SS"]])) {
    object[["SS"]]["Sigma"]
  } else cat("No data to extract variance-covariance matrix from", fill = TRUE)
}

summary.IntData <- function(object, digits = getOption("digits"),
  verbose = FALSE, ...) {

  cat("Variables in model: ", colnames(object$ModelData), "\n", fill = TRUE)
  cat("Test of Simple Slopes: \n\n")
  print(coef(object), digits = digits)
  cat("\nResidual DF: ", object[["SS"]][["df.residual"]], fill = TRUE)
  if (verbose) {
    cat("\n\n ~~~~ Verbose Output ~~~~")
    cat("\nRaw coefficients: \n")
    print(object[["SS"]]["B"], digits = digits)
    cat("\nVariance-Covariance Matrix: \n")
    print(vcov(object), digits = digits)
    cat("\nContrast Matrix: \n")
    print(object[["SS"]]["model"], digits = digits)
  }
}

SimpleSlopes <- function(object, x, intvar, intval, values) {
  ## Simple Slopes
  values[, x] <- 1

  tt <- terms(object)
  Terms <- delete.response(tt)
  intval <- list(unique(intval))
  names(intval) <- intvar
  tmp <- expand.grid(intval)
  slopedat <- do.call(rbind, rep(list(values), nrow(tmp)))
  slopedat[, colnames(tmp)] <- tmp

  MF <- model.frame(Terms, slopedat, na.action = "na.omit",
    xlev = object$xlevels)
  if (!is.null(cl <- attr(Terms, "dataClasses")))
    .checkMFClasses(cl, MF)
  X <- model.matrix(Terms, MF, contrasts.arg = object$contrasts)
  L <- attr(Terms, "factors")
  mode(L) <- "logical"
  zero.index <- !attr(X, "assign") %in% which(L[x, ])
  X[, zero.index] <- 0
  B <- matrix(coef(object))
  Sb <- vcov(object)
  slopes <- as.vector(X %*% B)
  SE <- sqrt(diag(X %*% Sb %*% t(X)))
  df <- object$df.residual
  tval <- slopes/SE
  pval <- 2 * pt(q = abs(tval), df = df, lower.tail = FALSE)
  coeff <- data.frame(slopes, SE, tval, pval)
  colnames(coeff) <- c("Simple Slope", "Std. Error", "t value", "Pr(>|t|)")

  output <- list(coefficients = coeff, B = B, Sigma = Sb, model = X, df.residual = df)
  return(output)
}

#tmp <- do.call("rbind", with(tout$PredictedData, by(cbind(mpg, hp), wt, FUN = mean)))
#tmp <- cbind(tmp, tout$coefficients)
#p +geom_text(data = tmp, aes(x = hp, y = ifelse(mpg < 0, mpg * .95, mpg * 1.05), label = format(`Simple Slope`, digits = 2)))
